E15test3

Date: 08.04.21

We use WGCNA to run an iterative analysis on a data set. If no features are prided, HVGs are calculated in the script in the script, using Seurat, the mvp method a threshold of 0.25 for the dispersion.

Pre-analysis

We need to first set up our working environment.

my.date = format(Sys.Date(), "%d.%m.%y")

# The following setting is important, do not omit
backup.options = options()
options(stringsAsFactors = FALSE)

Parse all the parameters and data we need to be able to run the analysis. We need the data in seurat format. This must contain at least 20 samples (no problems with sc data), according to the documentation of the package itself.

# placeholder for the variables and options we will get from an upper level pipeline

project_name = params$project.name

p.Wdata = params$sc.data

gnames = params$gene.names
gnames[,2]= as.character(gnames[,2])
rownames(gnames) = gnames[,1]

Wdata = params$data

# The subset in form of cell identities, F if using the whole sample
my.subset = params$cells
# Variable genes, or genes to use, F if the genes will be calculated in the script
my.vargenes = params$features

my.dr = params$reduction

dir.create(paste0(params$dir,"./WGCNA_results"), showWarnings = F)

my.sp = params$sp

is.pseudocell = params$is.pseudocell

my.min.cell = params$min.cell

Now we calculate the genes we will use.

if (params$GO==T) {
  
  if (my.vargenes == F){
    my.gouni = rownames(p.Wdata@assays$RNA@counts)[which(apply(p.Wdata@assays$RNA@counts, 1, function(x) length(which(x >0))) > my.min.cell-1)]
  } else {
      my.gouni = rownames(p.Wdata@assays$RNA@counts)[which(apply(p.Wdata@assays$RNA@counts, 1, function(x) length(which(x >0))) > 0)]
  }
}

# If we need to subset the data
if (my.subset) {
  p.Wdata = subset(p.Wdata, idents=my.subset)
}

# If no variable genes are provided
nonex = which(apply(p.Wdata@assays$RNA@counts, 1, function(x) length(which(x >0))) < my.min.cell)

if (my.vargenes == F) {
  # First get rid of non-expressed genes
  p.Wdata = subset(p.Wdata, features = rownames(p.Wdata@assays$RNA)[-nonex])
  
  # Find the variable genes
  p.Wdata=Seurat::FindVariableFeatures(p.Wdata, dispersion.cutoff = c(0.25,Inf), mean.cutoff = c(0,Inf), selection.method = "mvp", assay = "RNA")
  Seurat::VariableFeaturePlot(p.Wdata, assay = "RNA")
  Expr = Seurat::VariableFeatures(object = p.Wdata, assay = "RNA")
  
} else { Expr = my.vargenes }

# if (length(my.ortho) > 1) {
#   Expr = Expr[which(Expr %in% my.ortho[,1])]
#   Expr = Expr[ which(Expr %in% rownames(p.Wdata@assays$RNA)[-nonex]) ]
# }
 
if (is.pseudocell) {

  datExpr=Wdata@assays$RNA@counts[Expr,]
  # datExpr = datExpr[which(rowSums(as.matrix(datExpr))>0),]
  # Expr = rownames(datExpr)

} else{datExpr=Wdata@assays$RNA@data[Expr,]}

# Check the length
print(paste0("We have ", length(Expr), " genes in the variable genes object"))
## [1] "We have 2967 genes in the variable genes object"
# Check the size and transform
dim(datExpr)
## [1] 2967  513
datExpr = t(as.matrix(datExpr))

Now we need to calculate the soft threshold power. First it calculates the similarity and then transforms this similarity to a weighted network. The scale-free topology is calculated for each of the powers.

We choose the smallest power for which the scale-free topology fit index reaches 0.90. If none of the powers reaches 0.90, we take the one with the maximum, as long as we have a number above 0.75. If none of them reaches at least 0.75 we need to check our dataset.

## Warning: executing %dopar% sequentially: no parallel backend registered
## Warning in (function (x, y = NULL, robustX = TRUE, robustY = TRUE, use =
## "all.obs", : bicor: zero MAD in variable 'x'. Pearson correlation was used
## for individual columns with zero (or missing) MAD.
## Warning in (function (x, y = NULL, robustX = TRUE, robustY = TRUE, use =
## "all.obs", : bicor: zero MAD in variable 'y'. Pearson correlation was used
## for individual columns with zero (or missing) MAD.
##    Power SFT.R.sq   slope truncated.R.sq  mean.k. median.k. max.k.
## 1      1 0.000237  -0.608          0.610 1500.000  1.50e+03 1600.0
## 2      2 0.255000 -12.000          0.948  769.000  7.67e+02  895.0
## 3      3 0.147000  -4.970          0.812  400.000  3.95e+02  513.0
## 4      4 0.371000  -4.380          0.503  212.000  2.05e+02  310.0
## 5      5 0.913000  -5.200          0.889  115.000  1.09e+02  210.0
## 6      6 0.953000  -4.190          0.944   63.800  5.83e+01  154.0
## 7      7 0.969000  -3.380          0.965   36.700  3.16e+01  120.0
## 8      8 0.972000  -2.800          0.968   22.100  1.73e+01   98.0
## 9      9 0.972000  -2.400          0.964   13.900  9.64e+00   83.0
## 10    10 0.958000  -2.150          0.949    9.190  5.45e+00   72.0
## 11    12 0.946000  -1.800          0.946    4.670  1.82e+00   56.8
## 12    14 0.953000  -1.590          0.967    2.810  6.35e-01   46.4
## 13    16 0.952000  -1.470          0.970    1.890  2.31e-01   38.7
## 14    18 0.964000  -1.420          0.984    1.360  8.98e-02   32.8
## 15    20 0.956000  -1.400          0.976    1.030  3.64e-02   28.1
## 16    22 0.938000  -1.380          0.947    0.800  1.54e-02   24.2
## 17    24 0.921000  -1.370          0.924    0.637  6.57e-03   21.1
## 18    26 0.946000  -1.340          0.961    0.516  2.89e-03   18.4
## 19    28 0.962000  -1.320          0.985    0.425  1.31e-03   16.2
## 20    30 0.907000  -1.340          0.930    0.354  6.09e-04   14.3
##  [1]  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20

## [1] "Or power is 5"

WGCNA analysis

Now, the following piece of code will run WGCNA iteratively, to end up with out final modules. The iterations follow these steps:

  • Tree cutting for module calculation
    • Calculate an adjacency matrx from the data, then turn it into topological overlap and then into a distance matrix
    • Calculate the tree based on the topological overlap distance
    • Calculate the automatic height to cut out the 0.05 quantile
    • Generate a matrix where we calculate the amount of modules, size of de modules and wheter if we have a grey module based on:
      • Different minimum module sizes, arbitrarely set to 7:30
      • Different cut-heiughts going 0.0005 up and down from the automatic height in steps of 0.0001
    • Check if any combination of the parameters will get rid of the grey module
      • If we only get grey modules, subset the matrix for the height at which the grey module is the smallest
      • If we have a combination without grey module AND we have at least the same number of modules as in the beginning, we subset for whichever those heights are
    • Take wichever min module sizes gives us at least the same amount of clusters as in the beginning, if none, then the highest
    • Chose the max of the reminding min module sizes.
  • Calculate the actual modules
  • Calculate the eigengenes
  • Calculate the module membership per gene
  • Calculate the p.value of the membership to a given module
  • Get rid of the genes in the grey module
  • Delete the genes that are not significantely associated with their module
  • Save the remaining geneset
  • Print how many genes were deleted due to significance
  • Unless 0 genes were deleted in the last step, update the expression matrices and beginn again.

ONLY THE FIRST TWO ITERATIONS ARE DIFFERENT.

FIRST:

  • During the tree cutting
    • Set the minimum module size to arbitrary 15
    • Subset for the heights that get rid of at least 50% of the genes with that min module size
    • Choose the height that gives us the most clusters

SECOND:

  • Set the resulting number of modules as the ground number of modules
# Running WGCNA iteratively

my.Clnumber = 20
change = 0
genesets=list()
nonsig = 1

while(nonsig != 0) {
  
  my.adjacency =WGCNA::adjacency(datExpr,type = "signed", power = my.power, corFnc = "bicor")
  
  # Turn adjacency into topological overlap (high overlap if they share the same "neighborhood")
  TOM=WGCNA::TOMsimilarityFromExpr(datExpr,networkType = "signed", TOMType = "signed", power = my.power, corType = "bicor")
  
  #Put the names in the tree
  colnames(TOM) <- gnames[colnames(datExpr),2]
  
  rownames(TOM) <- gnames[colnames(datExpr),2]
  
  #Make it a distance
  dissTOM = 1-TOM
  
  # Call the hierarchical clustering function
  geneTree = hclust(as.dist(dissTOM), method = "average")
  
  
  # Here I calculate the cutting height. Using the same formula and approach that the WGCNA package uses for the automatic function
  nMerge = length(geneTree$height) # The whole height of the tree
  refQuantile = 0.05 # What's the quantile that we want to exclude
  refMerge = round(nMerge * refQuantile) 
  refHeight = geneTree$height[refMerge]
  cutheight = signif(0.99 * (max(geneTree$height) - refHeight) + refHeight,4)
  
  # We construct THE TABLE that will help us make decisions
  # Min cluster sizes, from 7 to 30
  x=seq(7,30,1)
  # The height, up and down from the calculated height. We expect some "No module detected"
  y=seq(cutheight-0.0005,cutheight + 0.0005,0.0001)
  
  # The actual dataframe
  w=data.frame()
  # Populate, with i=min cluster size, j=cutting height, z=total number of clusters, z.1.=what's the first cluster? 0 is grey 1 is something else,
  # z.1.'=what's the size of the first cluster?
  for (i in x) {
    for (j in y) {
      sink("aux")
      z=table(dynamicTreeCut::cutreeDynamic(dendro = geneTree,  method="tree", minClusterSize = i, deepSplit = T, cutHeight = j, verbose = 0))
      sink(NULL)
      v=data.frame(i,j,dim(z),names(z[1]),unname(z[1]))
      w=rbind(w,v)
    }
  }
  
  # The height is then the one where we have the least number of genes in the first cluster
  my.height = w$j[which(w$unname.z.1..==min(w$unname.z.1..))]
  
  # Since different heights can give us the minimum grey size, we chose the computed height, if present, or the highest one.
  if (cutheight %in% my.height) {
    my.Clsize = w[which(w$j == cutheight),]
  } else { my.Clsize = w[which(w$j == max(my.height)),] }
  
  
  # This is to know, if we're looking for a minimum of cluster numbers
  # If we still have a lot of genes, we don't want to limit the number of clusters
  # if ( ((dim(datExpr)[2]) / length(Expr)) > 0.6 ) {change = 0}
  
  # If this is the first iteration after 0.6 of the genes are gone, we assign the number of clusters (and an extra for the grey in the case)
  if (change == 2) {
    my.Clnumber = length(table(dynamicColors)) + 1
  }
  # Count another iteration
  change = change + 1
  
  # If we don't have a gray cluster anymore, then we subset for those rows, and set a new height. ONLY if we get the same amount of clusters!
  if (any(w$names.z.1.. == 1)) { #any combination gives us no grey
    
    my.Clsize = w[which(w$names.z.1.. == 1),] # Take all combinations that gives us no grey
    
    if (any(my.Clsize$dim.z. >= (my.Clnumber -1) )) { # If there is any giving us the determined amount or more
      my.Clsize = my.Clsize[which(my.Clsize$dim.z. >= (my.Clnumber - 1) ),,drop=F] # Subset for those
    } else { my.Clsize = my.Clsize[which(my.Clsize$dim.z. == max(my.Clsize$dim.z.)),,drop=F] } # Or for the highest
    
    # Take the ones with the smallest number of clusters
    my.Clsize = my.Clsize[which(my.Clsize$dim.z. == min(my.Clsize$dim.z.)),,drop=F]
    # Take the one with the highest min cluster size
    my.Clsize = my.Clsize[which(my.Clsize$i == max(my.Clsize$i)),,drop=F]
    
    if (cutheight %in% my.Clsize$j) { # if original computed height is in,
      my.height = cutheight # take it
    } else { my.height = max(my.Clsize$j) } # Otherwise, the highest height
    
    my.Clsize = max(my.Clsize$i)
    
  }
  
    # Subset the table again, for those sizes that will gives the same number of clusters or more. IF NONE, use the highest number
  if (!any(w$names.z.1.. == 1)){
    if (any(my.Clsize$dim.z. >= my.Clnumber)) {
      my.Clsize = my.Clsize[which(my.Clsize$dim.z. >= my.Clnumber),,drop=F]
    } else {
      my.Clsize = my.Clsize[which(my.Clsize$dim.z. == max(my.Clsize$dim.z.)),,drop=F]}
    
   # Take the ones with the smallest number of clusters
    my.Clsize = my.Clsize[which(my.Clsize$dim.z. == min(my.Clsize$dim.z.)),,drop=F]
    # Take the one with the highest min cluster size
    my.Clsize = my.Clsize[which(my.Clsize$i == max(my.Clsize$i)),,drop=F]
    
    if (cutheight %in% my.Clsize$j) { # if original computed height is in,
      my.height = cutheight # take it
    } else { my.height = max(my.Clsize$j) } # Otherwise, the highest height
    
    my.Clsize = max(my.Clsize$i)
  }
  
  # If we still have more than 60% of the genes, we just use the min size of 15 regardles
  if ( change < 3 ) {
    my.Clsize = 15
    my.height = cutheight
  }
  
  dynamicMods = dynamicTreeCut::cutreeDynamic(dendro = geneTree,  method="tree", minClusterSize = my.Clsize, deepSplit = T, cutHeight = my.height)
  
  #dynamicMods = cutreeDynamic(dendro = geneTree, distM = dissTOM, deepSplit = 2, minClusterSize = minModuleSize)
  
  table(dynamicMods)
  
  # Convert numeric lables into colors
  dynamicColors = WGCNA::labels2colors(dynamicMods)
  table(dynamicColors)
  
  # Plot the dendrogram and colors underneath
  WGCNA::plotDendroAndColors(geneTree, dynamicColors, "Modules",
                      dendroLabels = NULL,
                      cex.dendroLabels = 0.6,
                      addGuide = TRUE,
                      main = "Gene dendrogram and module colors",
                      guideAll = F)
  
  par(mfrow=c(1,1))
  
  # Calculate eigengenes
  MEList = WGCNA::moduleEigengenes(as.matrix(datExpr), colors = dynamicColors)
  
  MEs = MEList$eigengenes
  
  # Calculate the module membership
  geneModuleMembership = as.data.frame(WGCNA::signedKME(datExpr, MEs))
  MMPvalue = as.data.frame(WGCNA::corPvalueStudent(as.matrix(geneModuleMembership), nrow(datExpr)))
  
  # We're gonna make a list, where we keep the genes that are significantly associated with each module
  x=c()
  xy=list()
  
  # We also need a vector with all the dynamic colors
  dcols = 1:length(levels(as.factor(dynamicColors)))
  
  # Getting rid of the grey module
  grey.genes = length(which(dynamicColors == "grey"))
  if (any(levels(as.factor(dynamicColors)) == "grey")) {
    dcols = dcols[-which(levels(as.factor(dynamicColors)) == "grey")]
  }
  
  # Run the loop to get the genes
  for (i in dcols) {
    modGenes = rownames(MMPvalue)[which(dynamicColors==levels(as.factor(dynamicColors))[i] & MMPvalue[,i]<0.01)]
    x=c(x,modGenes)
    xy[[i]]=modGenes
    #print(paste0(levels(as.factor(dynamicColors))[i]," ",length(modGenes),
    #" of ", length(which(dynamicColors==levels(as.factor(dynamicColors))[i]))))
    #print(gnames[modGenes,2])
  }
  
  # Make a new list, where we keep ALL the gens thar are left from the iteration, that will be used to make the new object. To keep track
  genesets[[length(genesets)+1]] = colnames(datExpr)
  
  # Give me a message saying how many genes are gone this time
  cat( paste0( grey.genes, " genes not assigned to any module.", '\n',
                 length(which(!(colnames(datExpr)%in%x))) - grey.genes, " genes excluded due to significance."))
  # Save this also, cause if it's 0 then we stop the whole thing
  nonsig = length(which(!(colnames(datExpr)%in%x)))
  
  # If it ain't 0, subset the dynamic colors and the expression data
  if (length(which(!(colnames(datExpr)%in%x))) != 0) {
    dynamicColors=dynamicColors[-which(!(colnames(datExpr)%in%x))]
    datExpr=datExpr[,-(which(!(colnames(datExpr)%in%x)))]
  } 

}
## Warning in (function (x, y = NULL, robustX = TRUE, robustY = TRUE, use =
## "all.obs", : bicor: zero MAD in variable 'x'. Pearson correlation was used
## for individual columns with zero (or missing) MAD.
## TOM calculation: adjacency..
## ..will not use multithreading.
##  Fraction of slow calculations: 0.000000
## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.
## 0 genes not assigned to any module.
## 528 genes excluded due to significance.
## Warning in (function (x, y = NULL, robustX = TRUE, robustY = TRUE, use =
## "all.obs", : bicor: zero MAD in variable 'x'. Pearson correlation was used
## for individual columns with zero (or missing) MAD.

## TOM calculation: adjacency..
## ..will not use multithreading.
##  Fraction of slow calculations: 0.000000
## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.
## 0 genes not assigned to any module.
## 135 genes excluded due to significance.
## Warning in (function (x, y = NULL, robustX = TRUE, robustY = TRUE, use =
## "all.obs", : bicor: zero MAD in variable 'x'. Pearson correlation was used
## for individual columns with zero (or missing) MAD.

## TOM calculation: adjacency..
## ..will not use multithreading.
##  Fraction of slow calculations: 0.000000
## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.
## 0 genes not assigned to any module.
## 161 genes excluded due to significance.
## Warning in (function (x, y = NULL, robustX = TRUE, robustY = TRUE, use =
## "all.obs", : bicor: zero MAD in variable 'x'. Pearson correlation was used
## for individual columns with zero (or missing) MAD.

## TOM calculation: adjacency..
## ..will not use multithreading.
##  Fraction of slow calculations: 0.000000
## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.
## 0 genes not assigned to any module.
## 37 genes excluded due to significance.
## Warning in (function (x, y = NULL, robustX = TRUE, robustY = TRUE, use =
## "all.obs", : bicor: zero MAD in variable 'x'. Pearson correlation was used
## for individual columns with zero (or missing) MAD.

## TOM calculation: adjacency..
## ..will not use multithreading.
##  Fraction of slow calculations: 0.000000
## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.
## 0 genes not assigned to any module.
## 36 genes excluded due to significance.
## Warning in (function (x, y = NULL, robustX = TRUE, robustY = TRUE, use =
## "all.obs", : bicor: zero MAD in variable 'x'. Pearson correlation was used
## for individual columns with zero (or missing) MAD.

## TOM calculation: adjacency..
## ..will not use multithreading.
##  Fraction of slow calculations: 0.000000
## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.
## 0 genes not assigned to any module.
## 9 genes excluded due to significance.
## Warning in (function (x, y = NULL, robustX = TRUE, robustY = TRUE, use =
## "all.obs", : bicor: zero MAD in variable 'x'. Pearson correlation was used
## for individual columns with zero (or missing) MAD.

## TOM calculation: adjacency..
## ..will not use multithreading.
##  Fraction of slow calculations: 0.000000
## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.
## 0 genes not assigned to any module.
## 8 genes excluded due to significance.
## Warning in (function (x, y = NULL, robustX = TRUE, robustY = TRUE, use =
## "all.obs", : bicor: zero MAD in variable 'x'. Pearson correlation was used
## for individual columns with zero (or missing) MAD.

## TOM calculation: adjacency..
## ..will not use multithreading.
##  Fraction of slow calculations: 0.000000
## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.
## 0 genes not assigned to any module.
## 3 genes excluded due to significance.
## Warning in (function (x, y = NULL, robustX = TRUE, robustY = TRUE, use =
## "all.obs", : bicor: zero MAD in variable 'x'. Pearson correlation was used
## for individual columns with zero (or missing) MAD.

## TOM calculation: adjacency..
## ..will not use multithreading.
##  Fraction of slow calculations: 0.000000
## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.
## 0 genes not assigned to any module.
## 2 genes excluded due to significance.
## Warning in (function (x, y = NULL, robustX = TRUE, robustY = TRUE, use =
## "all.obs", : bicor: zero MAD in variable 'x'. Pearson correlation was used
## for individual columns with zero (or missing) MAD.

## TOM calculation: adjacency..
## ..will not use multithreading.
##  Fraction of slow calculations: 0.000000
## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.
## 0 genes not assigned to any module.
## 1 genes excluded due to significance.
## Warning in (function (x, y = NULL, robustX = TRUE, robustY = TRUE, use =
## "all.obs", : bicor: zero MAD in variable 'x'. Pearson correlation was used
## for individual columns with zero (or missing) MAD.

## TOM calculation: adjacency..
## ..will not use multithreading.
##  Fraction of slow calculations: 0.000000
## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.

## 0 genes not assigned to any module.
## 0 genes excluded due to significance.

We calculate here the expression of genes using the single-cell data

p.MEList = MEList

raw.datExpr = p.Wdata@assays$RNA@data[colnames(datExpr),]

raw.datExpr = t(as.matrix(raw.datExpr))

raw.MEList = WGCNA::moduleEigengenes(raw.datExpr, colors = dynamicColors)
  
p.MEList = raw.MEList

Modules of co-expression

For single cells We can see what are the expression levels of our co-expression modulesin the single cells. Using the dimensionality reduction provided. Default is a tSNE

 xx=list()
 yy=c(levels(as.factor(dynamicColors)))
for (i in 1:length(yy)) {
  toplot = data.frame(Seurat::Embeddings(p.Wdata[[my.dr]]))

  xx[[i]] = ggplot2::ggplot(toplot[order(p.MEList$averageExpr[,i]),],
                   ggplot2::aes_string(x=colnames(toplot)[1], y=colnames(toplot)[2])) +
    ggplot2::geom_point(ggplot2::aes_string(color=p.MEList$averageExpr[order(p.MEList$averageExpr[,i]),i]), size=2) +
    ggplot2::scale_size(range = c(1, 1)) +
    ggplot2::theme_void() +
    ggplot2::theme(legend.position="none") +
    ggplot2::scale_colour_gradientn(colours = c("gray90", "gray90", yy[i], yy[i])) +
    ggplot2::labs(colour=levels(as.factor(dynamicColors))[i])
}
gridExtra::grid.arrange(grobs=xx, ncol=4)

Create a list object with all the essential data from the analyses

WGCNA_data = list()
WGCNA_data[["datExpr"]] = datExpr
WGCNA_data[["dynamicMods"]] = dynamicMods
WGCNA_data[["MEList"]] = MEList
WGCNA_data[["MEs"]] = MEs
WGCNA_data[["modGenes"]] = xy
WGCNA_data[["MEList.sc"]] = p.MEList
WGCNA_data[["genesets"]]= genesets
WGCNA_data[["TOM"]]= TOM
WGCNA_data[["adjacency"]]= my.adjacency

saveRDS(WGCNA_data, file = paste0(params$dir,project_name,"_WGCNA_data_",my.date,".rds"))

To create the visualizations of the actual co-expression networks, we use this code. It also creates files that can be read into cytoscape

# we create a new directory where the network files will be
my.dirname = paste0(params$dir,project_name,"WGCNA_networks_",my.date,"/")

dir.create(my.dirname)
## Warning in dir.create(my.dirname): '/scicore/home/tschoppp/feregrin/test/
## scWGCNA_testing/E15test3WGCNA_networks_08.04.21' already exists
modules=unique(dynamicColors) # The modules

# To generate the files necessary for the network visualization in cytoscape
for (i in 1:length(modules)) {
  mod=is.finite(match(dynamicColors, modules[i])) # Each module
  cyt = WGCNA::exportNetworkToCytoscape(TOM[mod,mod],
                                 edgeFile = paste( paste0(my.dirname,"CytoscapeInput-edges-"),
                                                  paste(modules[i], collapse="-"), ".txt", sep=""), # The name of the edge file
                                 nodeFile = paste( paste0(my.dirname,"CytoscapeInput-nodes-"),
                                                  paste(modules[i], collapse="-"), ".txt", sep=""), # The node file
                                 weighted = TRUE, #Puts the weights
                                 threshold = 0.0, #Threshold for adjacency, 0=all genes
                                 nodeNames = colnames(datExpr)[mod], #The names of the nodes/genes
                                 altNodeNames = gnames[colnames(datExpr)[mod],2], #Other names
                                 nodeAttr = geneModuleMembership[colnames(datExpr)[mod],i]) #Some more info about the nodes
}

# We build a list where we keep the network plots
netplots=list()

lcols = c(rep("black", length(unique(dynamicColors))))
lcols[(apply(col2rgb(levels(as.factor(dynamicColors))), 2,
               function(x) (x[1]*0.299 + x[2]*0.587 + x[3]*0.114)) < 75)] = "white"

mynets = list()

# A loop that makes the network plots in R
for (i in 1:length(xx)) {
  # Load directly edges tables from the file we created
  mynetwork = read.table(paste0(my.dirname,"CytoscapeInput-edges-",levels(as.factor(dynamicColors))[i],".txt"),
                       header = T, stringsAsFactors = F, fill=T)
  # We get rid of all the edges that have a very small weight. The cutoff set to keep at least ONE edge on the nodes
  x = max( c( min(aggregate(mynetwork$weight, by = list(mynetwork$fromNode), max)$x),
              min(aggregate(mynetwork$weight, by = list(mynetwork$toNode), max)$x) ) )
  mynetwork=mynetwork[-which(mynetwork$weight < x),]

  # We rescale the weights so that we have them from 0 to 1
  mynetwork$weight01 = GGally::rescale01(mynetwork[,3])
  # And then multiply them by 2, to give the heavy edges a 2 thickness, we ad 0.2 to not have completelly invisible edges
  mynetwork[,3] = (GGally::rescale01(mynetwork[,3]) * 2) + 0.2
  # Convert this edgelist into a network object
  mynet = network::network(mynetwork, matrix.type="edgelist", ignore.eval=F)

  # Load in the nodes table
  mynodes = read.table(paste0(my.dirname,"CytoscapeInput-nodes-",levels(as.factor(dynamicColors))[i],".txt"),
                     header = T, stringsAsFactors = F, fill = T)
  rownames(mynodes) <- mynodes$nodeName
  # put in the membership value, from the ModuleMembership table, and rescale that from 0 to 1.
  mynodes$membership = GGally::rescale01(geneModuleMembership[rownames(mynodes),i])
  # We multiply by 30 to give a good range of sizes, plus 1 to avoid innexistant nodes
  mynodes$membership= (mynodes$membership*30)+1

  mynet = network::set.vertex.attribute(mynet, "membership", mynodes[network::network.vertex.names(mynet),"membership"])

  mynets[[i]] = mynet

}

To check the results module by module, we report some GO terms and individual plots, with module sizes and gene names. First the GO terms analyses

if (params$GO==T) {
  
  #Get the ENSEMBL to ENTREZ list from the BioMart

  require(paste0("org.",my.sp,".eg.db"),character.only = T)
  my.ENS2EG = get(paste0("org.",my.sp,".egENSEMBL2EG"), pos = paste0("package:org.",my.sp,".eg.db"))
  
  IDs=as.list(my.ENS2EG)

  #Create two lists, one for the DEGs and one for the GOs
  dewg = list()
  gowg = list()

  # A loop that goes through the DEG files we created earlier
  for(i in 1:length(yy)){ #as many clusters as we have
     dewg[[i]] = xy[[i]] #Read the table
     ex = sapply(dewg[[i]], function(x) AnnotationDbi::exists(x, my.ENS2EG)) #Which genes have an ENTREZ?
    dewg[[i]] = dewg[[i]][ex] #Only those genes
    dewg[[i]] = unlist(IDs[dewg[[i]]]) #The ENTREZ IDs
    gowg[[i]] = limma::goana(dewg[[i]], species = my.sp, universe = unlist(IDs[my.gouni])) #The GO analysis
   }
  
  
  
  goplots=list()
  goterms=list()

  for (i in 1:length(yy)) {
    toplot=limma::topGO(gowg[[i]], n= 50, ontology = c("BP"))

    #Change to character and then back to factor, to keep the order from TopGO
    toplot$Term = as.character(toplot$Term)
    toplot$Term = factor(toplot$Term, levels = unique(toplot$Term))

    goterms[[i]]=toplot

  }

  #TO make the acutal plot we're having, we combine all the GOs
  g.plot = data.frame()
  topterms = data.frame()

  for (i in 1:length(goterms)) {
    topterms = goterms[[i]][1:5,]
    topterms$cluster = i
    g.plot = rbind(g.plot, topterms)
  }

  g.plot$cluster = as.factor(g.plot$cluster)

  g.plot$Term = droplevels(g.plot$Term)
  my.terms = nchar(levels(g.plot$Term))>45
  my.gterms = strtrim(levels(g.plot$Term), 45)
  my.gterms[my.terms] = paste0(my.gterms[my.terms],"...")
  levels(g.plot$Term) = my.gterms
  
  
  my.terms = nchar(as.character(g.plot$Term))>45
  my.gterms = strtrim(g.plot$Term, 45)
  my.gterms[my.terms] = paste0(my.gterms[my.terms],"...")

  ggplot2::ggplot(g.plot, ggplot2::aes(x=Term, y=-log10(P.DE), fill=cluster)) +
      ggplot2::geom_point(ggplot2::aes(shape = cluster), size = 3) +
      ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90, hjust = 1)) +
      ggplot2::labs(x="GO Term", y="-log10 p. value") +
      ggplot2::scale_fill_manual(name = "Module",
                        labels = yy,
                        values = yy) +
      ggplot2::scale_shape_manual(name = "Module",
                         labels = yy,
                         values = rep(c(21,22,23,24,25),20)) +
    ggplot2::coord_flip()
}

And here the report of each of the modules

# Go trhough the modules
for(i in 1:length(xx)){

  # A table is done, where we put the genes names that are making-up the module, in rows of 10. We fill the las row with empty spaces
  ktable1 = data.frame( matrix( c(gnames[xy[[i]],2][order(geneModuleMembership[xy[[i]],i], decreasing = T)],
                                  rep(" ", 10-length(gnames[xy[[i]],2])%%10 )), ncol = 10, byrow = T ) )
  colnames(ktable1) = rep(".", 10)

  if (params$GO == T) {
    # Another table, where we put the top 10 GO terms, in two columns
    ktable2 = data.frame( goterms[[i]][1:10,])
    ktable2[,5] = as.character(signif(ktable2[,5],4))
    colnames(ktable2) = c("Term", "Ont", "N", "n", "Adj. p-value")
    }

  # We make a plot in which we show the mean expression level of the co-expression module

  toplot = data.frame(Seurat::Embeddings(p.Wdata[[my.dr]]))

  myplot = ggplot2::ggplot(toplot[order(p.MEList$averageExpr[,i]),],
                   ggplot2::aes_string(x=colnames(toplot)[1], y=colnames(toplot)[2])) +
      ggplot2::geom_point(ggplot2::aes_string(color=p.MEList$averageExpr[order(p.MEList$averageExpr[,i]),i])) +
      ggplot2::scale_size(range = c(1, 1)) +
      ggplot2::theme_void() +
      ggplot2::scale_colour_gradientn(colours = c("gray90", "gray90", yy[i], yy[i])) +
      ggplot2::labs(colour=levels(as.factor(dynamicColors))[i])

  mynet=mynets[[i]]

  mynetplot = GGally::ggnet2(mynet,
                         mode = "fruchtermanreingold",
                         layout.par = list(repulse.rad=network::network.size(mynet)^1.1,
                                           area=network::network.size(mynet)^2.3), # Give space in the middle
                         node.size = network::get.vertex.attribute(mynet,'membership'), max_size = 20, #The size of the nodes
                         node.color = levels(as.factor(dynamicColors))[i], # The color of the module
                         edge.size = "weight",
                         edge.color = "black",
                         edge.alpha = network::get.edge.attribute(mynet,'weight01'),
                         ) +
    ggplot2::theme(legend.position="none") +
    ggplot2::geom_label(ggplot2::aes(label=gnames[network::network.vertex.names(mynet),2]),
               fill = levels(as.factor(dynamicColors))[i],
               alpha = 0.5,
               color=lcols[i],
               fontface = "bold")

  # We put our mean expression plot toghether with the network plot
  gridExtra::grid.arrange(grobs=list(myplot,mynetplot), nrow=2)
  # # plot.new()
  # # dev.off()
  cat("\n")
  # A table showing the number, color and size of the module
  print(knitr::kable(data.frame(module=i, color=names(table(dynamicColors))[i], size=unname(table(dynamicColors)[i]))))
  # cat("\n")
  # The other two tables we just made
  print(knitr::kable(ktable1))
  if (params$GO == T) {
   print(knitr::kable(ktable2))
  }
  cat("\n")


}

module color size
1 black 105
. . . . . . . . . .
Krt5 Perp Krt15 Sfn Krt14 Trim29 Fgfbp1 Lgals7 Wnt6 Trp63
Mt2 Pkp3 Bcam Epcam Dsp Wnt7b Gjb6 Esrp1 Aqp3 Gjb2
Kremen2 Col17a1 Dmkn Pdlim1 Rab38 Lmo1 Fermt1 Pkp1 Serpinb10 Wnt3
Wnt10a Gata3 Meis2 Mapk13 Tfap2a Ckmt1 Tfap2b Cdh1 Ptprf Lypd6b
Serpinb5 Irf6 Cstdc5 Anxa8 Tfap2c Krt17 Barx2 Krtdap Fzd10 Irx4
Igsf9 Hspb8 Lamc2 Irx2 Zfp296 Fras1 Itga3 Klc3 Ddr1 Gm13219
Cebpa Esrp2 Frem2 Ppp1r14c Fzd6 Endou Chchd10 Lamb3 Prxl2a Hotairm1
Wnt4 Car12 Ifi202b Pltp Them5 St14 Hoxb4 Etl4 Cgn Cldn6
Egfr Jag1 Adh1 Plet1 Ces2g Mgst1 Nrarp Urah Ch25h Trim2
Dpp4 Bicdl2 Akap17b Atp8b1 Pard6b Prss27 Rdh12 Tubb3 Nipal1 Spns2
Gchfr Anxa9 Tmem184a Gprc5a Ly6g6d
Term Ont N n Adj. p-value
GO:0060429 epithelium development BP 910 32 3.103e-14
GO:0048729 tissue morphogenesis BP 570 26 3.575e-14
GO:0002009 morphogenesis of an epithelium BP 472 24 3.857e-14
GO:0009888 tissue development BP 1556 38 5.893e-12
GO:0043588 skin development BP 209 15 3.082e-11
GO:0008544 epidermis development BP 229 15 1.128e-10
GO:0009887 animal organ morphogenesis BP 864 26 4.227e-10
GO:0048513 animal organ development BP 2634 44 1.823e-08
GO:0009913 epidermal cell differentiation BP 132 10 4.322e-08
GO:0007155 cell adhesion BP 1010 25 5.033e-08
## Warning: Removed 3 rows containing missing values (geom_label).

module color size
2 blue 226
. . . . . . . . . .
Crabp1 Egfl6 Tmem132c Prrx2 Vcam1 Mmp11 Twist1 Map1b Efna5 Tril
Nkd2 Prrx1 C1qtnf2 Irx3 Dpysl3 Col7a1 Ngfr a Adamts2 Zeb2
Col6a3 Fez1 Lmna Rspo1 Ctsk Septin11 Serpinc1 Irx1 Cntfr Ednra
Sp5 Tmem119 Apela Col23a1 Frem1 Lix1 Col5a1 Ccn4 Col14a1 Ptk7
Grem2 Eva1b Cpxm2 Adam33 Tbx18 Rasl11b Plxna4 Emx2 Sema5a Creb3l1
Fam174b Calb1 Nt5e Cbr3 Ltbp1 Cd248 Hmgn3 Cpne5 Pdgfrb Kcnk2
Cdo1 Cxcl14 Zfp503 Ppp1r14a Vegfd Itih5 Tac1 Gpc4 Ephb2 Spats2l
Nradd Rgs7bp Dnm3 Thy1 Vstm4 Wif1 Gm14964 St8sia2 Fndc1 Lrig1
C1qtnf7 En1 C130021I20Rik Alx3 Axin2 Stmn2 Itga8 Egflam Thsd7a Hic1
Pax1 Olfml3 Zfp536 Lama2 Cdc42ep5 Fgf10 Igfbp2 Bend5 Rhoj Rnd3
Slit3 Igsf10 Dkk1 Unc5d Wnt11 Cox20 Svbp Hunk Slc7a10 Mme
Lrrn3 Kcna4 Cnih2 Tpbg Lin7a Tspan11 Septin7 Sdc2 Pappa2 Dpt
Chst2 Serping1 Scube2 Tgfbr3 Kcns3 Ism1 Frmd6 Col8a1 Mex3b Mylip
Vstm2b Nnmt Tex15 Ptprd Cspg5 Iglon5 Rai14 Cttnbp2 Smarca2 Hand2
Cacna2d3 Bmp3 Amot Palm Plpp3 Cacng4 Rgs17 Adam22 Cyth3 Plppr3
Fnbp1l Dtx4 Bex4 Scara3 Hand2os1 Lmx1b Homer2 NA Cplx2 Inhba
Samd14 Rbp1 Lsamp Fjx1 Pdlim4 Fap Gsta4 Pcdh11x Cald1 Frmd4a
Aspn Robo1 Rdh10 Vcan Adamts5 Wbp1 Stra6 Stox2 Dchs1 NA
Enox1 Fgf18 Scx Bex1 Fbn1 Ptchd4 Grik5 Uchl1 Sema3a Cdh2
Insyn1 Fat4 Shd H2bu2 Fblim1 Klhdc8b Tmem26 Tnmd Rarres2 Postn
Pdlim3 Ebf2 Rimklb Gstt1 Sbk1 Mkx Adgrg2 Tpm2 Epha3 Fn1
Thbs2 Kirrel3 Chodl Gulp1 Pianp Pdpn Gm10371 Pnck NA Tia1
Prr16 Ccdc136 Sbspon Tmeff2 Pdgfc Tmem30b
Term Ont N n Adj. p-value
GO:0009653 anatomical structure morphogenesis BP 2239 94 5.912e-22
GO:0032501 multicellular organismal process BP 5039 147 1.41e-21
GO:0048731 system development BP 3637 122 2.079e-21
GO:0007275 multicellular organism development BP 4034 127 4.338e-20
GO:0048856 anatomical structure development BP 4376 132 1.711e-19
GO:0032502 developmental process BP 4644 134 4.243e-18
GO:0009887 animal organ morphogenesis BP 864 51 3.86e-17
GO:0072359 circulatory system development BP 971 53 2.329e-16
GO:0051239 regulation of multicellular organismal process BP 2519 87 1.223e-14
GO:0048869 cellular developmental process BP 3314 102 3.195e-14
## Warning: Removed 1 rows containing missing values (geom_label).

module color size
3 brown 223
. . . . . . . . . .
Cd52 Coro1a Psmb8 Bin2 Laptm5 Rac2 Hcls1 Cd53 Evi2a Myo1g
Apbb1ip Fermt3 Selplg Vav1 H2-D1 Ptprc Mndal Napsa Srgn Plac8
Lcp1 Psmb9 Wfdc17 Fmnl1 Nrros Fxyd5 Itgb2 Ccl6 Gpsm3 Nckap1l
Vsir Cytip Fcgr2b Fgd2 Lat2 Arhgap9 Gmfg Sp140 Cyp4f18 Ncf4
Ighm Cfp H2-K1 Cd68 Casp1 Ptpro Lcp2 AI662270 Sh2d1b2 Slc15a3
C3ar1 Cx3cr1 Ifitm6 Ptpn18 Cd44 Cxcl2 Lrmp Tnfaip8l2 Otulinl Cd74
Sp110 Clec4d Bcl2a1b Lpxn Gm2a Pira2 Ccl3 Rgs2 Fyb Dok2
Mcemp1 Ifi203 Pkib Ppfia4 Ifi211 Ly9 Plbd1 Siglece C1qb Fes
Dock2 Emb Il1b Tm6sf1 C1qa Hfe Tpd52 Cd14 Lyn C5ar1
Ltb4r1 Rgs10 Trf Arhgap15 Hck Tifab Tlr2 Rgs18 Capg Cxcl3
Glipr1 Cd33 S100a4 Osgin1 Slc11a1 Tifa Sirpa Zc3hav1 C1qc Syngr2
Cd79b Rtp4 Rgs1 Klhl6 Ccr5 Adcy7 Oasl2 Phf11b Irf8 Rab7b
Ramp1 Plaur Ptpn7 Cd300lf Gm14548 Il4ra Arhgap45 Skap2 Bcl2a1d Fcrls
Mcub NA H2-Aa Lgals9 Bank1 H2-DMa Clec5a Ccl4 Pirb Cysltr1
Rassf4 Milr1 H2-Eb1 Slamf9 Atp2a3 Ccl2 Ptprj Clec10a Il13ra1 Nlrp3
Pmaip1 Ehd4 Csf3r Ms4a7 Hacd4 Fgl2 Tmem37 Clec2i Tuba4a B4galnt1
Tmem106a Hpgds Adssl1 Cd83 H2-Oa Il1r2 Wfdc21 Emilin2 Il18rap Gpr65
Slc9a3r1 Bst2 Fam129a Bmp2k Gch1 Mir22hg Prtn3 Rab32 Ahr Stom
Xlr Hhex H2-Ab1 Apobec3 Selenop Pglyrp1 Adap2os Adgre5 Fam107b Ccrl2
F13a1 Lpcat2 Gpr137b-ps Dusp5 Osm Prkcb Tnni2 Klrd1 P2ry12 5430427O19Rik
Gas7 Pak1 H2-T23 Bcl2a1a Blnk Hmox1 Mrc1 Ang Trpv2 Hp
Cstdc4 S100a8 Pf4 Smim3 Hsd11b1 Atp1a3 H2-DMb1 S100a9 Mgl2 Gbp7
Aph1c Kcnn4 F10
Term Ont N n Adj. p-value
GO:0006955 immune response BP 908 98 1.446e-58
GO:0002376 immune system process BP 1714 122 6.129e-55
GO:0006952 defense response BP 966 81 7.363e-39
GO:0002682 regulation of immune system process BP 997 75 1.739e-32
GO:0002684 positive regulation of immune system process BP 683 63 5.928e-32
GO:0050776 regulation of immune response BP 540 56 6.835e-31
GO:0045087 innate immune response BP 467 52 3.901e-30
GO:0045321 leukocyte activation BP 644 59 1.163e-29
GO:0001775 cell activation BP 729 62 2.041e-29
GO:0050778 positive regulation of immune response BP 429 47 7.127e-27

module color size
4 cyan 58
. . . . . . . . . .
Ms4a6c Ctss Pld4 Clec4a3 Mpeg1 Cybb Ly86 Clec4a2 Igsf6 Lyz2
Clec4a1 Ms4a6b Adgre1 Ebi3 AB124611 Trem2 Alox5ap Apoc2 Slfn2 Clec12a
Cd300c2 Ms4a6d Clec4n Cd48 Lst1 Fcgr3 Ccr2 Ctsc Tnfrsf13b Ms4a4c
Plek Ifi27l2a Unc93b1 Csf1r Ifi204 Ncf2 Ccr1 Apobec1 Cyth4 Ccl9
Tyrobp Lrrc25 Irf5 Spi1 Cd86 Gm49339 Fcgr1 Ptafr Fcer1g P2ry6
Ifi207 Rnase6 Epsti1 Clec7a Pilra Ptpn6 Myo1f Aif1
Term Ont N n Adj. p-value
GO:0006952 defense response BP 966 28 4.625e-18
GO:0006955 immune response BP 908 27 1.244e-17
GO:0002376 immune system process BP 1714 34 2.37e-17
GO:0009605 response to external stimulus BP 1557 28 9.315e-13
GO:0045087 innate immune response BP 467 17 2.309e-12
GO:0051707 response to other organism BP 590 18 8.772e-12
GO:0043207 response to external biotic stimulus BP 592 18 9.28e-12
GO:0009607 response to biotic stimulus BP 617 18 1.844e-11
GO:0002274 myeloid leukocyte activation BP 158 11 3.032e-11
GO:0006954 inflammatory response BP 494 16 6.654e-11

module color size
5 green 151
. . . . . . . . . .
Top2a Birc5 Cdca8 Nusap1 Ccna2 Hmgb2 Cdca3 Ube2c Tpx2 Smc4
Tacc3 Spc25 Aurkb Spc24 Racgap1 Prc1 Cks2 Pbk Lockd Tuba1b
Smc2 H2az1 Kif23 Mki67 Hmmr Pclaf Ccnb1 Incenp Ube2s Kif22
Cenpe Stmn1 Cdk1 Plk1 Cenpf Nuf2 Cdc20 Kifc1 Tubb5 Tubb4b
H2ax Mis18bp1 H2az2 Melk Ccnb2 Bub3 Ndc80 Cenpa Bub1b Knstrn
Ncapd2 Depdc1a Pimreg Lmnb1 Tk1 Cdkn2c Arl6ip1 Rrm2 Aurka Kif20a
Ckap2l Asf1b H1f10 Shcbp1 Bub1 Sapcd2 Ccdc34 Nucks1 Kif20b Psrc1
Hmgn2 Mxd3 Sgo2a Neil3 Dlgap5 Anln Tyms Knl1 Ncapg Tmpo
Ckap2 Nde1 Kif2c Kpna2 Cdkn3 Hist1h2ap Cit Cdc25c Dek Esco2
Calm2 Jpt1 Ube2t Ptma Tubb6 Gmnn Aspm Rad51 H2ac8 H4c9
Tuba1c Dut Hmgb3 Plk4 Pkmyt1 Gas2l3 Atad2 H1f0 Lig1 Cbx3
Tcf19 H1f4 Uhrf1 Cep89 Syce2 Lsm6 Gm47283 Fxyd6 Rpgrip1 Mcm5
Supt16 Slbp 1500009L16Rik Gins2 Slc24a5 Pcna Chaf1b Wdhd1 Xist Gli1
Ptms Slc29a1 Hells Hs6st2 Mid1 Pik3r3 Syne2 Sms Hoxc10 Clasp2
Nsg1 Spock3 Mcm3 Ccne1 Lnpk Hoxa9 Fhl1 Fzd4 Mcm2 Tafa4
Col4a2
Term Ont N n Adj. p-value
GO:0007049 cell cycle BP 1408 92 2.361e-54
GO:0000278 mitotic cell cycle BP 730 67 1.762e-46
GO:0022402 cell cycle process BP 972 74 3.145e-46
GO:1903047 mitotic cell cycle process BP 596 60 1.872e-43
GO:0051301 cell division BP 516 56 4.151e-42
GO:0000280 nuclear division BP 339 42 1.779e-33
GO:0048285 organelle fission BP 380 42 2.117e-31
GO:0140014 mitotic nuclear division BP 235 35 1.23e-30
GO:0051276 chromosome organization BP 956 58 1.056e-29
GO:0007059 chromosome segregation BP 275 36 1.817e-29

module color size
6 greenyellow 70
. . . . . . . . . .
Tspo Gpx1 Nfkbia Rap1b Msn Actr3 Cdk2ap2 Ier5 Lrrfip1 Prr29
Serpina3g Lat Il22 Sec11c Tgfb1 Cd4 Tox2 Zap70 Ifngr1 Ikzf2
Cited4 Nrgn Arg1 Il17re Cyfip2 Septin1 Cited2 Hpcal1 Tec Ifi30
Actg1 Macroh2a1 Dusp1 Espn Arl4c Tspan13 Gsn Samhd1 Jmjd1c Trim8
Cmtm7 Pik3r1 Cdkn1a Ifrd1 Zfp36 Gsto1 Samd10 Hexb Mzb1 Cnn2
Fam124b Sdf2l1 Car2 Ctsd Zyx Ahnak Akr1c13 Rrad Anxa1 Dnase2a
Myb Kdm7a Ptgs2 Cxcl16 Ly6a Slc7a11 Prr11 Crip1 Asb17 Serp1
Term Ont N n Adj. p-value
GO:0002376 immune system process BP 1714 32 5.504e-12
GO:0002682 regulation of immune system process BP 997 21 1.029e-08
GO:0006955 immune response BP 908 19 6.893e-08
GO:0034097 response to cytokine BP 637 16 7.865e-08
GO:0070887 cellular response to chemical stimulus BP 2093 29 9.315e-08
GO:0006950 response to stress BP 2695 33 1.464e-07
GO:0006952 defense response BP 966 19 1.823e-07
GO:0071310 cellular response to organic substance BP 1686 25 3.017e-07
GO:0010033 response to organic substance BP 2075 28 3.075e-07
GO:0050896 response to stimulus BP 5705 50 3.383e-07

module color size
7 grey60 40
. . . . . . . . . .
Mfap4 Fstl1 Mdk Gas1 Mfap2 Col3a1 Col1a1 Fbn2 Cdh11 Pdgfra
Tcf4 Il11ra1 Mmp14 6330403K07Rik Fbln5 Macroh2a2 Col1a2 Igfbp4 Sfrp2 Col6a2
Gas2 Tnfaip6 Runx1t1 Clmp Wnt5a Cxcl12 Bgn Loxl1 Nfib Ubb
Fgfr1 Hoxd12 Sox11 Pfn2 Marcksl1 Zfp36l1 Tmsb10 Vim Gm9844 Tmsb4x
Term Ont N n Adj. p-value
GO:0001501 skeletal system development BP 465 14 8.579e-12
GO:0009888 tissue development BP 1556 20 4.213e-10
GO:0009653 anatomical structure morphogenesis BP 2239 23 7.685e-10
GO:0060322 head development BP 553 13 1.194e-09
GO:0030334 regulation of cell migration BP 768 14 6.247e-09
GO:0030198 extracellular matrix organization BP 229 9 8.173e-09
GO:0035295 tube development BP 948 15 1.03e-08
GO:2000145 regulation of cell motility BP 800 14 1.052e-08
GO:0007417 central nervous system development BP 692 13 1.778e-08
GO:0040012 regulation of locomotion BP 854 14 2.412e-08

module color size
8 lightcyan 50
. . . . . . . . . .
Cyp11a1 Il4 Cpa3 Il1rl1 Cd200r3 Cma1 Hdc Prss34 Alox5 Skap1
Mcpt8 Cd200r4 Adora3 Tpsg1 Tpsb2 Hs3st1 Cst7 Gata2 Cela1 Smpx
Tph1 Slc45a3 Rab27b Tpsab1 Il13 Kit Gnaz Man2b1 Prkcq Gadd45a
Tesc Gyg B2m Slc7a5 Ftl1 Plin2 Cebpb Cyba Klk8 Ubash3b
Nr4a2 Ly6e Pim1 Snx2 Tiparp Metrnl Fth1 Psap Sat1 Gstt2
Term Ont N n Adj. p-value
GO:0006954 inflammatory response BP 494 15 9.618e-11
GO:0032101 regulation of response to external stimulus BP 601 13 1.208e-07
GO:0006952 defense response BP 966 16 1.238e-07
GO:0043303 mast cell degranulation BP 37 5 1.875e-07
GO:0002279 mast cell activation involved in immune response BP 38 5 2.154e-07
GO:0002448 mast cell mediated immunity BP 38 5 2.154e-07
GO:0033008 positive regulation of mast cell activation involved in immune response BP 17 4 3.284e-07
GO:0043306 positive regulation of mast cell degranulation BP 17 4 3.284e-07
GO:0150076 neuroinflammatory response BP 47 5 6.425e-07
GO:0033005 positive regulation of mast cell activation BP 21 4 8.171e-07

module color size
9 lightgreen 31
. . . . . . . . . .
Bcl11b Tnfrsf19 Megf6 Ass1 Ank3 Fzd10os S100a11 Acaa2 Fst Gja1
Tom1l1 Ccdc3 Sult5a1 Cdh3 Wnt16 Ltbp4 Anxa2 Efcab15 1110028F18Rik Susd4
Svil Mpped1 Dlx4os Calm4 Ociad2 Sox15 Dlx2 Fam83f Cldn1 Limk2
Sbsn
Term Ont N n Adj. p-value
GO:0043616 keratinocyte proliferation BP 38 4 1.26e-06
GO:0043588 skin development BP 209 6 4.888e-06
GO:0008544 epidermis development BP 229 5 0.0001212
GO:0021877 forebrain neuron fate commitment BP 9 2 0.0001618
GO:0009888 tissue development BP 1556 11 0.0002589
GO:0010644 cell communication by electrical coupling BP 12 2 0.0002954
GO:0030855 epithelial cell differentiation BP 448 6 0.0003435
GO:0042475 odontogenesis of dentin-containing tooth BP 67 3 0.0004008
GO:0060429 epithelium development BP 910 8 0.0005477
GO:2000810 regulation of bicellular tight junction assembly BP 17 2 0.0006046

module color size
10 lightyellow 22
. . . . . . . . . .
Mlana Tyr Dct Slc45a2 Pax3 Pmel Ptgds Gpnmb Foxd3 Mcoln3
Ednrb Trpm1 Plp1 Syngr1 Insc Bace2 Sox10 Npy Gsta1 Nkain4
Phyhipl Slc26a7
Term Ont N n Adj. p-value
GO:0043473 pigmentation BP 91 7 7.477e-11
GO:0048066 developmental pigmentation BP 48 5 1.099e-08
GO:0042438 melanin biosynthetic process BP 19 4 1.93e-08
GO:0006582 melanin metabolic process BP 21 4 2.974e-08
GO:0044550 secondary metabolite biosynthetic process BP 24 4 5.263e-08
GO:0046189 phenol-containing compound biosynthetic process BP 37 4 3.226e-07
GO:0019748 secondary metabolic process BP 46 4 7.895e-07
GO:0046148 pigment biosynthetic process BP 47 4 8.62e-07
GO:0042440 pigment metabolic process BP 56 4 1.758e-06
GO:0018958 phenol-containing compound metabolic process BP 69 4 4.081e-06
## Warning: Removed 1 rows containing missing values (geom_label).

module color size
11 magenta 85
. . . . . . . . . .
S100a14 Cldn3 Sprr1a Gm16136 Ppl Mab21l4 Rab25 Nectin4 Tacstd2 Aldh3b2
Bspry Cldn4 Lypd3 Upk1b Krt4 Krt19 Asprv1 Ovol1 Areg Paqr5
Krt6a Cldn23 Mal Sult2b1 Dsc2 Clic3 Gpr87 Krt7 Fam25c Ly6d
Upk3bl Grhl3 Spint1 Ly6g6e Rhov Krt18 Adtrp Cers3 Hspb1 Grhl1
Pdzk1ip1 Crb3 Sptlc3 Ephx3 Sp6 Dsg1a Phldb3 Map3k6 Evpl Tmem125
Krt8 Cnfn Ablim1 Slc2a3 Capns2 Gng13 Nebl Smtnl2 Tmem40 Jup
Cst6 Mall Csta1 NA Scel Tent5b Il1rn Paqr6 Zfp750 Snx31
Kctd1 Mpzl3 Eps8l2 Lor Sytl1 Psca Aldh1a7 Prss8 Arap2 Krt1
Gipc2 Tmem79 Dsc3 Klf5 Prom2
Term Ont N n Adj. p-value
GO:0043588 skin development BP 209 17 3.483e-15
GO:0008544 epidermis development BP 229 16 2.644e-13
GO:0031424 keratinization BP 15 7 1.148e-12
GO:0009913 epidermal cell differentiation BP 132 12 1.46e-11
GO:0030216 keratinocyte differentiation BP 82 10 4.228e-11
GO:0030855 epithelial cell differentiation BP 448 18 8.352e-11
GO:0060429 epithelium development BP 910 23 1.379e-09
GO:0061436 establishment of skin barrier BP 17 5 3.607e-08
GO:0033561 regulation of water loss via skin BP 22 5 1.5e-07
GO:0009888 tissue development BP 1556 26 4.686e-07

module color size
12 midnightblue 54
. . . . . . . . . .
mt-Cytb mt-Atp6 mt-Co3 mt-Co2 mt-Nd1 Fabp5 mt-Co1 mt-Nd2 Pard6g Sostdc1
Tmem132a Npnt Spry1 1700093K21Rik mt-Nd5 Sema3f Cachd1 Cygb Ddit4 Hspa1a
Nr2e3 Cyb5a Lypd6 Irx6 Dcn Rhbdl3 Hoxc8 Ier3 Astn2 Bmp7
Msx2 Edar Jag2 Mfhas1 Emp2 Syt1 Hoxc13 Ptprz1 Gdpd1 Foxo6
Prr5l Nav2 Notch3 Dkk4 Unc5b Tnfrsf12a F2rl1 Syt13 Narf Cldn7
Arhgdib Kctd11 5730522E02Rik Zfp185
Term Ont N n Adj. p-value
GO:0046034 ATP metabolic process BP 211 9 6.985e-08
GO:0042773 ATP synthesis coupled electron transport BP 60 6 8.257e-08
GO:0009205 purine ribonucleoside triphosphate metabolic process BP 234 9 1.691e-07
GO:0009126 purine nucleoside monophosphate metabolic process BP 236 9 1.818e-07
GO:0009167 purine ribonucleoside monophosphate metabolic process BP 236 9 1.818e-07
GO:0009199 ribonucleoside triphosphate metabolic process BP 238 9 1.953e-07
GO:0009161 ribonucleoside monophosphate metabolic process BP 240 9 2.097e-07
GO:0009144 purine nucleoside triphosphate metabolic process BP 242 9 2.25e-07
GO:0009123 nucleoside monophosphate metabolic process BP 247 9 2.675e-07
GO:0009141 nucleoside triphosphate metabolic process BP 257 9 3.74e-07

module color size
13 pink 89
. . . . . . . . . .
Cd37 Trbc2 Hcst Ptprcap Itgal Cd69 Ltb Il7r Il18r1 Cd7
Klrb1b Il2rg Sla Nkg7 Itgb7 Ccl5 Gpr183 Cd160 Ikzf1 Parvg
Itk Rasal3 Trbc1 Cd3g Tnfaip3 Tnfrsf25 Tagap Tnfrsf9 Cxcr6 Tcf7
Samsn1 Fam83a Gpr171 Rora Xlr4b Xlr4c Slc6a13 Xlr4a Ctsw Lta
Rab37 Tmem154 Icos Cd28 Asb2 Serpinb1a Rhoh Pfn1 Sh3bgrl3 Crem
Npl Icam1 Cd47 Tagln2 Ckb Ccr7 Flt3 Klf6 Esyt1 Itm2b
Ecm1 Ifitm1 Rassf2 Arpc1b Irak2 Abca3 Serinc3 Slc16a10 Clic1 Actb
Junb Cxcl10 Csrnp1 Arl6ip5 Dleu2 Rbm38 Slc25a5 S100a13 Rnf130 Akr1b3
Rnasel Fam111a Zfp263 Ctnnd2 Apoc1 Runx1 Entpd1 Pou2f2 Mitf
Term Ont N n Adj. p-value
GO:0002376 immune system process BP 1714 44 1.186e-18
GO:0042110 T cell activation BP 373 21 4.866e-15
GO:0002682 regulation of immune system process BP 997 31 7.02e-15
GO:0046649 lymphocyte activation BP 538 24 7.118e-15
GO:0001775 cell activation BP 729 27 9.683e-15
GO:0006955 immune response BP 908 29 3.485e-14
GO:0045321 leukocyte activation BP 644 25 4.162e-14
GO:0002521 leukocyte differentiation BP 432 20 9.331e-13
GO:0002252 immune effector process BP 526 21 3.933e-12
GO:0046631 alpha-beta T cell activation BP 122 12 8.744e-12

module color size
14 purple 72
. . . . . . . . . .
Hoxd13 Ccnd2 Sdc1 Foxp2 Ntng1 Ccnd1 Creb5 Hoxa11os Gdf5 Serinc2
Adgrl3 Hoxd11 Dkk3 Fbln2 Mpped2 Col12a1 Pgrmc1 Aopep Epha4 Pthlh
Hhip Msx1 Car14 H3f3b Foxp4 Tob1 Prdm16 Sema3c Hmcn1 Lhx2
Tgfb3 Ptn Igsf3 Serpine2 Col8a2 Serpinf1 Erg Mn1 Lmo4 Hes1
Tgfb2 Fxyd7 Ssbp2 Syt11 Thsd4 Mecomos Pbx1 Hoxa10 Nrxn1 Nrtn
Frk Kcnd2 Axl Bmp1 Fam171b Erc2 Nell2 Plxdc2 Tenm3 Car11
Kcnmb4 Dlg5 Srsf12 Sv2a Tmtc2 2700069I18Rik Ank2 Mllt3 Htra1 Nrn1
Slit2 Pcdh10
Term Ont N n Adj. p-value
GO:0048513 animal organ development BP 2634 39 2.857e-12
GO:0009888 tissue development BP 1556 28 2.396e-10
GO:0048732 gland development BP 356 14 1.466e-09
GO:0035295 tube development BP 948 21 2.388e-09
GO:0032502 developmental process BP 4644 47 2.662e-09
GO:0009887 animal organ morphogenesis BP 864 20 3.003e-09
GO:0032501 multicellular organismal process BP 5039 49 3.003e-09
GO:0030154 cell differentiation BP 3150 38 3.688e-09
GO:2000026 regulation of multicellular organismal development BP 1758 28 3.964e-09
GO:0050793 regulation of developmental process BP 2188 31 6.499e-09

module color size
15 red 119
. . . . . . . . . .
Actc1 Tnnt2 Tnnt1 Cdh15 Mymx Myog Chrna1 Mymk Vgll2 Des
Chrng Rbm24 Neb Atp2a1 Ank1 Arpp21 Ttn Myod1 Fitm1 Cryab
Tnni1 Myl1 Pnmal2 Gatm Jsrp1 Rtn2 Megf10 Mylpf Acta2 Klhl41
Hspb2 Tnik Cap2 Nes Iffo1 Tex14 Fndc5 Mrln Myl4 Traf3ip3
Il17b Pitx2 Pgam2 Srpk3 Mstn Srl Rgs16 Olfml2b Hdac11 Myh3
Sln Flnc Car3 Tnnc2 Cck Tnnc1 Celf2 Myl6b Tubb2b Septin4
Ralgps2 Igfbp3 Prkaca Pdlim7 Pgm5 Mycl Acta1 Casq2 Dll1 Ssc5d
Parm1 Tppp3 Cox6a2 Unc45b Casz1 Ccdc155 Tmem35a Nexn Fmo1 Arhgdig
Lrrn1 Kcne1l Tceal6 Eya2 Cited1 AW551984 Gpr37 Dbndd1 Six1 Frmd4b
Tln2 Msc Vash2 Pdgfa Tmem200a Rxrg Shisal2a Hspd1 Ptgis Fam110b
Spon1 Itga7 Trim63 Myf5 Shisa2 Rai2 Lurap1l Epor Tcea3 Mad2l2
Sypl2 Necab1 Bcl7a Rhobtb1 Dbx1 Rplp1 Dcx Cbln2 Camk2n1
Term Ont N n Adj. p-value
GO:0061061 muscle structure development BP 585 38 4.37e-24
GO:0007517 muscle organ development BP 350 28 4.561e-20
GO:0003012 muscle system process BP 312 26 4.263e-19
GO:0014706 striated muscle tissue development BP 373 27 3.078e-18
GO:0060537 muscle tissue development BP 394 27 1.246e-17
GO:0006936 muscle contraction BP 217 20 1.113e-15
GO:0060538 skeletal muscle organ development BP 170 18 2.901e-15
GO:0006941 striated muscle contraction BP 123 16 4.153e-15
GO:0007519 skeletal muscle tissue development BP 163 17 2.335e-14
GO:0051146 striated muscle cell differentiation BP 264 20 4.833e-14

module color size
16 salmon 65
. . . . . . . . . .
Crabp2 Lgals1 Ddah2 Marcks Lpar1 Tle5 Nbl1 Col26a1 Col5a2 Basp1
Rbms3 Lhfp Epb41l3 Nnat Fbln1 Tead2 Bcl11a Dbn1 Podxl2 Tpm1
Col6a1 S1pr3 Chd3 Tuba1a Phldb2 Ctxn1 Ntf3 Cacnb3 Glipr2 Dclk1
Pcdh18 Myh10 Ogn Dkk2 Nfia Lhfpl2 Ly6h Foxp1 Csrp2 Lox
Myo1b Zfhx4 Sncaip Angptl4 Sulf1 Sox4 Cped1 Adamtsl1 Efnb3 Osr2
Robo2 Osr1 Mir100hg Bmp4 Cacna1g Aplp1 Nfatc4 Csrp1 Agtr2 Crlf1
Cmtm3 Epha7 Mecom Igfbp5 Col4a1
Term Ont N n Adj. p-value
GO:0048856 anatomical structure development BP 4376 46 2.725e-11
GO:0007275 multicellular organism development BP 4034 44 3.949e-11
GO:0009653 anatomical structure morphogenesis BP 2239 33 4.221e-11
GO:0072001 renal system development BP 251 13 7.913e-11
GO:0009888 tissue development BP 1556 27 1.798e-10
GO:0032502 developmental process BP 4644 46 2.579e-10
GO:0001655 urogenital system development BP 289 13 4.53e-10
GO:0032501 multicellular organismal process BP 5039 47 1.119e-09
GO:0048731 system development BP 3637 39 3.535e-09
GO:0007399 nervous system development BP 1857 26 4.442e-08
## Warning: Removed 1 rows containing missing values (geom_label).

module color size
17 tan 66
. . . . . . . . . .
Ddit4l Mmp13 Gpx3 Rbp4 Cp Ffar4 Ibsp Steap4 Fabp3 Slc36a2
Enpp6 Ocm AW112010 Itgb3 Fam20c Cemip Maf Agt Hpgd Bmp8a
2200002D01Rik Dapk2 Phex Mgst2 Serpini1 Hand1 Tcim Lpl Frzb Mmp9
Loxl4 Chst1 Lratd1 Plaat3 Sp7 Irx5 Gm10638 Ptges Clu Cdkn2b
Csf1 Hmgcs2 Tnfaip2 Entpd3 Cds1 B3galt2 A230083N12Rik Podnl1 9130024F11Rik Dmp1
Syna Dlx5 Cd55 Foxq1 Rcan2 Slpi Pla2g4a Lgr6 Gm19705 Rasl11a
Jam2 NA Ackr3 Nos1ap Adamts9 Ank
Term Ont N n Adj. p-value
GO:0001503 ossification BP 330 11 2.063e-07
GO:0009888 tissue development BP 1556 21 2.042e-06
GO:0043062 extracellular structure organization BP 263 9 2.366e-06
GO:1900020 positive regulation of protein kinase C activity BP 7 3 2.944e-06
GO:1900019 regulation of protein kinase C activity BP 7 3 2.944e-06
GO:0030198 extracellular matrix organization BP 229 8 7.758e-06
GO:0061448 connective tissue development BP 251 8 1.513e-05
GO:0031214 biomineral tissue development BP 126 6 2.035e-05
GO:0048513 animal organ development BP 2634 26 2.805e-05
GO:0010743 regulation of macrophage derived foam cell differentiation BP 16 3 4.576e-05

module color size
18 turquoise 308
. . . . . . . . . .
Col9a3 Col9a1 Col2a1 Col11a1 Col9a2 Mia Hapln1 Cnmd Matn1 Acan
Col27a1 Col11a2 Matn3 Fkbp9 Susd5 Cthrc1 Plod2 Rcn3 Snorc Matn4
Papss2 Chadl P4ha1 Comp Serpinh1 Fibin Tspan4 Sox9 Wwp2 Scrg1
Dlk1 Itm2a Ugdh S100b Gale Kdelr3 Loxl3 Slc16a3 Fkbp10 Prelp
Islr Vkorc1 Fgfr3 Rflna Rcn2 Peg3 C1qtnf3 Chad Cspg4 H19
Plagl1 Fbln7 Higd1a Pantr1 Mgp Bnip3 Rtl3 Sox5 Snhg18 Pcolce2
Meltf Kdelr2 Cdkn1c Thbs1 Cytl1 Kazald1 P4ha2 Nfatc2 Cpe Maged2
Epyc Ccn1 P3h3 Cd9 Ier3ip1 Fzd9 Rian S100a1 Meg3 Neat1
Gpc3 Fos Sorbs2 Trim47 Fkbp11 Anxa5 Nrk Klf2 Sdc4 Mfge8
Capn6 Slc26a2 Pxylp1 Rlbp1 Cd24a Sfrp1 Il17d Kcnq1ot1 Gdf10 Tmem97
Scd2 Asb4 Pmp22 Dap Selenom Nid2 4930523C07Rik Scin Fosb Timp1
Fmod Phlda2 Sparc Mest Crispld1 Fam180a Aldoc Phyh Ndufa4l2 Jund
Dhx58os Lgals3 Lbhd2 Egr1 Srm Scd1 Jun Nog Cgref1 Zbtb20
Msmo1 Slc16a4 Fam162a Lmcd1 Cox4i2 Clec11a Ldhb Mif Ctsz Kcnk1
Alkbh1 Bex2 Steap1 Pitx1 Ppa1 Pcolce Tspan8 Klf4 Vwc2 Emilin1
Ndrg2 Moxd1 Tsc22d1 BC006965 Auts2 Serpina3n Ntm Tgfbi Ucma Tox
Tgm2 Hspa5 Lum Prkg2 Ebf1 Rgs3 Rgcc Foxc1 Fdps Tubb2a
Tbx15 Shox2 Dhcr7 Calml3 Klf9 Insig1 Timp3 Ier2 Hmgcr Pxdc1
Angptl1 Enpp1 Id4 Tspan32 Bcl2 Col15a1 Srpx Cpm Tent5a Fam89a
Prss35 Ooep Klhl13 Lpar4 1700086L19Rik Saa1 Pdcd4 Chil1 Arl4a Ecrg4
Fgfr2 Smoc2 Aebp1 Phlda1 Id3 Scn1b Hmga2 Gm26532 Hsp90b1 Uba7
Lrig3 Ihh Bhlhe41 Dbi Bambi Cyp51 Idi1 Bmp2 Pam Trabd2b
Rspo3 Cox6b2 Ppp2r2b Smim5 Ddah1 Spsb4 Dusp14 H1f2 Nim1k Klk10
Slc1a3 Lss Boc Cpxm1 Gem Bmper Saa2 2810403D21Rik Runx3 Smpd3
Fabp7 Mgarp Foxc2 Cdc42ep3 Hmgcs1 Panx3 Runx2 Tshz2 Bmp5 Car9
Acat2 Trib1 Smim14 Pth1r Sqle Enpp2 Aacs Mxra8 Clic4 2610528A11Rik
Efna1 Mvd Wfdc12 Tmem158 Zfp24 Barx1 Ung Tcf7l2 Serinc5 Htra3
Sgms2 Cpq Id1 Fam81a Ppic Raph1 Id2 Ttll3 Fbp1 Ccn2
Stk32b Smoc1 Arsi Mpz 1700049E15Rik Psma8 Cd82 Cdh10 Pcp4 Cyp26b1
Trps1 Ifitm5 Uts2b Batf3 3830403N18Rik Ltbp3 Ostn Angptl2 Efcab10 Tspan18
Flrt2 Btbd3 Penk Pcsk6 Sema3d Cldn10 Dab2 Gng11
Term Ont N n Adj. p-value
GO:0061448 connective tissue development BP 251 41 1.028e-24
GO:0051216 cartilage development BP 190 35 5.215e-23
GO:0001501 skeletal system development BP 465 51 2.159e-22
GO:0001503 ossification BP 330 43 7.132e-22
GO:0030198 extracellular matrix organization BP 229 35 3.146e-20
GO:0009888 tissue development BP 1556 90 1.978e-19
GO:0043062 extracellular structure organization BP 263 36 3.787e-19
GO:0002062 chondrocyte differentiation BP 114 25 1.19e-18
GO:0048705 skeletal system morphogenesis BP 235 31 3.709e-16
GO:0048731 system development BP 3637 142 8.386e-16

module color size
19 yellow 213
. . . . . . . . . .
Cdh5 Ecscr Icam2 Rasip1 Kdr Cldn5 Tie1 Pecam1 Myct1 Gimap6
Emcn Kank3 Ctla2a Plvap Plxnd1 Cd93 Sox17 Gpihbp1 Aplnr Tmem255a
Grrp1 Robo4 Esam Cav1 Adgrl4 Lmo2 Col18a1 Thsd1 Npr1 Cd34
Gimap4 Mfng Klhl4 Fam167b Eng Vamp5 Sox7 Rbpms Afap1l1 Flt1
Arap3 Bik Gngt2 Adgrf5 Myzap C130074G19Rik Hlx Apold1 Ipo11 Tm4sf1
Tek Fabp4 Dusp2 Smagp Prcp Arhgef15 Arhgap18 Cavin2 Clec1a Crmp1
Flt4 Dok4 Ushbp1 Dll4 Ptprb Sipa1 Mmrn2 Gimap5 Scarf1 Grap
Madcam1 Ly6c1 Abi3 Crip2 Bok S100a16 Dysf N4bp3 Gimap1 Podxl
Acvrl1 Bcl6b F11r Pcdh17 Gm45837 Trp53i11 Mmp15 Lpar6 Tfpi Tspan15
Tmem88 Cd36 Nova2 Apln Fgd5 Cav2 Gm16104 Sh3tc1 Cracr2b Cavin3
Ehd2 BC028528 Prkch Fkbp1a Sult1a1 Gap43 Dusp3 Stab1 Adam15 Cmtm8
Esm1 Cda C1qtnf9 Tpm4 Epas1 Ppp1r16b Mmrn1 Ets1 Pdgfb Ece1
Fam43a Ccm2l Tnfaip8l1 Adm Hdac7 Cd38 Cdh13 Rftn1 Tmem252 Kcne3
Gm13889 Tinagl1 Akap12 Npdc1 Ndrg1 Pcdh1 Procr Ccl21a Nectin2 Nrp2
Kitl Gja4 Cngb1 Cd40 Hid1 Nts Gm525 Palmd Clec14a Lamb1
Rasgrp2 Prrg3 Cavin1 Lck Kcna5 Ldb2 Calcrl Hoxb6 Mcam 4930542C12Rik
Itga6 Arhgap31 Prox1 Lyve1 Reln Lama4 Zfp979 Hbegf Depp1 Ephb1
Cd1d1 Arhgap29 Tcf15 Vwa1 Cotl1 Angpt2 F2r Slc16a13 Nedd9 Stc1
Nr2f2 Enpp3 AU021092 Tc2n Sema6d Gpm6a Hoxd9 Ccl21d Arhgef28 Fam171a2
Snrk Ifitm3 Meox2 Cpa1 Meox1 Filip1 Igf1 Hoxb9 Marchf11 Pgf
Fhl2 Edn1 Ccdc149 Apoa2 Myo1c Arhgef7 Cnrip1 Tspan12 Exoc3l4 Vwf
Tmod2 Lhx6 Adcy10
Term Ont N n Adj. p-value
GO:0001944 vasculature development BP 646 56 7.329e-28
GO:0072358 cardiovascular system development BP 657 56 1.738e-27
GO:0048514 blood vessel morphogenesis BP 538 51 4.094e-27
GO:0001568 blood vessel development BP 617 54 4.934e-27
GO:0001525 angiogenesis BP 445 47 5.207e-27
GO:0048646 anatomical structure formation involved in morphogenesis BP 966 60 2.903e-22
GO:0072359 circulatory system development BP 971 60 3.777e-22
GO:0035239 tube morphogenesis BP 786 53 3.528e-21
GO:2000145 regulation of cell motility BP 800 51 2.685e-19
GO:0040012 regulation of locomotion BP 854 52 8.254e-19
options(backup.options)

sessionInfo()
## R version 3.6.0 (2019-04-26)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: CentOS Linux 7 (Core)
## 
## Matrix products: default
## BLAS/LAPACK: /scicore/soft/apps/OpenBLAS/0.3.1-GCC-7.3.0-2.30/lib/libopenblas_sandybridgep-r0.3.1.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
## [1] org.Mm.eg.db_3.8.2   AnnotationDbi_1.46.1 IRanges_2.18.0      
## [4] S4Vectors_0.22.1     Biobase_2.44.0       BiocGenerics_0.30.0 
## [7] scWGCNA_0.1.0       
## 
## loaded via a namespace (and not attached):
##   [1] backports_1.1.4       Hmisc_4.2-0           sn_1.6-1             
##   [4] plyr_1.8.4            igraph_1.2.4.1        lazyeval_0.2.2       
##   [7] splines_3.6.0         listenv_0.7.0         ggplot2_3.1.1        
##  [10] TH.data_1.0-10        robust_0.4-18         digest_0.6.25        
##  [13] foreach_1.4.4         htmltools_0.3.6       GO.db_3.8.2          
##  [16] gdata_2.18.0          magrittr_1.5          checkmate_1.9.3      
##  [19] memoise_1.1.0         fit.models_0.5-14     cluster_2.0.9        
##  [22] doParallel_1.0.14     ROCR_1.0-7            limma_3.40.0         
##  [25] sna_2.5               globals_0.12.4        fastcluster_1.1.25   
##  [28] matrixStats_0.54.0    sandwich_2.5-1        colorspace_1.4-1     
##  [31] rrcov_1.4-7           blob_1.1.1            ggrepel_0.8.2        
##  [34] xfun_0.20             dplyr_1.0.2           crayon_1.3.4         
##  [37] jsonlite_1.6.1        impute_1.58.0         survival_2.44-1.1    
##  [40] zoo_1.8-5             iterators_1.0.10      ape_5.3              
##  [43] glue_1.4.0            gtable_0.3.0          leiden_0.3.3         
##  [46] future.apply_1.2.0    DEoptimR_1.0-8        scales_1.0.0         
##  [49] mvtnorm_1.0-10        GGally_1.4.0          DBI_1.0.0            
##  [52] bibtex_0.4.2          Rcpp_1.0.1            metap_1.3            
##  [55] plotrix_3.7-5         viridisLite_0.3.0     htmlTable_1.13.1     
##  [58] reticulate_1.15       foreign_0.8-71        bit_1.1-14           
##  [61] rsvd_1.0.3            preprocessCore_1.46.0 Formula_1.2-3        
##  [64] tsne_0.1-3            htmlwidgets_1.3       httr_1.4.0           
##  [67] gplots_3.0.1.1        RColorBrewer_1.1-2    acepack_1.4.1        
##  [70] TFisher_0.2.0         ellipsis_0.3.0        Seurat_3.1.4         
##  [73] ica_1.0-2             reshape_0.8.8         pkgconfig_2.0.2      
##  [76] nnet_7.3-12           uwot_0.1.8            dynamicTreeCut_1.63-1
##  [79] labeling_0.3          tidyselect_1.1.0      rlang_0.4.7          
##  [82] reshape2_1.4.3        munsell_0.5.0         tools_3.6.0          
##  [85] generics_0.0.2        RSQLite_2.1.1         statnet.common_4.3.0 
##  [88] ggridges_0.5.1        evaluate_0.13         stringr_1.4.0        
##  [91] yaml_2.2.0            npsurv_0.4-0          knitr_1.22           
##  [94] bit64_0.9-7           fitdistrplus_1.0-14   robustbase_0.93-5    
##  [97] caTools_1.17.1.2      purrr_0.3.4           RANN_2.6.1           
## [100] pbapply_1.4-1         future_1.13.0         nlme_3.1-140         
## [103] compiler_3.6.0        rstudioapi_0.11       plotly_4.9.2.1       
## [106] png_0.1-7             lsei_1.2-0            tibble_3.0.3         
## [109] pcaPP_1.9-73          stringi_1.4.3         highr_0.8            
## [112] lattice_0.20-38       Matrix_1.2-17         multtest_2.38.0      
## [115] vctrs_0.3.4           mutoss_0.1-12         pillar_1.4.6         
## [118] lifecycle_0.2.0       Rdpack_0.11-0         lmtest_0.9-37        
## [121] RcppAnnoy_0.0.16      data.table_1.12.2     cowplot_1.0.0        
## [124] bitops_1.0-6          irlba_2.3.3           gbRd_0.4-11          
## [127] patchwork_1.0.0       R6_2.4.0              latticeExtra_0.6-28  
## [130] network_1.16.0        KernSmooth_2.23-15    gridExtra_2.3        
## [133] codetools_0.2-16      MASS_7.3-51.4         gtools_3.8.1         
## [136] sctransform_0.2.1     mnormt_1.5-5          multcomp_1.4-10      
## [139] grid_3.6.0            rpart_4.1-15          coda_0.19-2          
## [142] tidyr_1.1.2           rmarkdown_2.6         Rtsne_0.15           
## [145] numDeriv_2016.8-1     WGCNA_1.67            base64enc_0.1-3